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The two-flavour Schwinger model is used to test the local boson action algorithm of M Liischer. The autocor- 
relation time is found to rise linearly with the number of auxiliary boson fields. An extension to the algorithm is 
proposed which exactly reproduces the full dynamical partition function for any number of boson fields. 



On 
O 

Q 



00 

o 
o 

<N 

OS 

Oh 

<D 



1. Introduction 

1.1. The Schwinger model 

The Schwinger model (2d QED) with stag- 
gered fermions provides a computationally simple 
testing ground for dynamical fermion algorithms. 
The theory is asymptotically free, so results relat- 
ing to algorithm performance should be similar to 
simulations of QCD. 

The model provides a useful probe of the low- 
eigenvalue behaviour of the algorithm which is 
particularly critical for the local boson method 
since the polynomial approximation breaks down 
here. The presence of approximate zero modes [3] 
(AZMs) associated with the topologically charged 
sectors of the theory prove difficult to handle. 

1.2. Liischer's local bosonic algorithm 

The algorithm [1] exploits a bosonic path in- 
tegral to simulate the fermion determinant of the 
full theory via a polynomial approximation to the 
inverse of the fermion matrix. The partition func- 
tion for the Schwinger model with 2 flavours of 
fermions is 



VU detQ e 



(1) 



Q is an hermitian version of the fermion matrix, 
scaled so it has eigenvalues in the range (—1, 1) 
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(2) 



Consider an n"^ order polynomial approximation, 
Vn{s) to 1/s where < s < 1. The polynomial 
has n roots in the complex plane, and can be 
written 



Vn{s) ocY[{s - Zi) 



(3) 



The optimal choice of roots comes from an anal- 
ysis of Chebyshev polynomials. The lower scale 
for accelerated convergence is set by the choice of 
parameter, e. 

With a set of roots, {z^} the determinant of Q 
can be replaced by 



detQ^ 



detP„(Q=) nr=i det(Q= - Zi) 



(4) 



For the two flavour theory, det P„((5^) is re- 
quired. The hermitian fermion matrix of (2) has 
eigenvalues in ± pairs which leads to a degener- 
acy in pairs of determinants in the product of (4). 
The square root is taken by removing one of the 
determinants in each pair. 

A set of auxiliary bosons is introduced with 
partition function, Zla which mimics the fermion 
determinant 



Zla[U] = j VxfVx)e- 



x'fiQ-y^)iQ-/^)xt 



(5) 



so that det Q Zla 

Note that the number of auxiliary fields is half 
the order of the polynomial for two flavour simu- 
lations. 

1.3. Implementation 

Updates are as follows: 

• the boson fields are updated from a (gaus- 
sian) heatbath, since they have a local 
quadratic action. 
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Figure 1. The plaquette integrated autocorrela- 
tion time vs. order of polynomial (e = 0.005) 



• gauge field updates are over-relaxation only. 
Ergodicity is ensured as the total action of 
the theory is updated stochastically when 
the boson fields are updated. There is lit- 
tle advantage in stochastic updates of the 
gauge fields. 



2. Autocorrelation times 

The integrated autocorrelation time for the pla- 
quette is shown as a function of n and e for a 
16 X 16 lattice, bare fermion mass = 0, /? = 3.0. 

From figure 1 we see that Tint oc f^- The 
CPU autocorrelation time (assuming the boson 
updates dominate) will therefore rise in propor- 
tion to . An estimate for the number of fields 
required to simulate fermions of mass m, is n oc 
1/m. This implies the cost of the algorithm will 
rise at least as 1/m?. 

Figure 2 shows that the autocorrelation time 
also rises as the lower scale of the approximation 
is reduced. The polynomial roots lie on an ellipse 
with focii at e and 1 and with minor axis length 
2y^. As e is reduced, the smallest eigenvalues of 
each boson coupling matrix will decrease like e so 
the autocorrelation time will rise approximately 
like 1/e 

The increase is not alleviated by increasing the 
denominator of (2). For staggered fermions, (un- 
like Wilson fermions [2]) the spectrum contains 
few eigenvalues near to s = 1. 
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Figure 2. The plaquette integrated autocorrela- 
tion time vs. e (n = 80) 

3. LARD - Local Action -|- Reduced De- 
terminant 

The full partition function of the theory is 
VU det Qe 



Employing the partition function of the local 
bosonic theory, (5), 



VU det Q JdetV{Q") Zla[U] e 



(7) 



If the polynomial approximation were exact, the 
product of the two determinants in (7) would be 
1. In practise, use of a computationally accessible 
number of polynomial terms (= number of boson 
fields) means there is some error in the polyno- 
mial and hence fluctuation in the product over 
configurations. The product, Ord quantifies the 
error in the polynomial on every configuration. V 
is a polynomial, hence it is simultaneously diago- 
nisable with Q so 



V/2 



(8) 



where are one of each of the degenerate pairs of 
eigenvalues oiQ The reduced determinant, Ord, 
has the following properties: 

1. the local boson method is "quenching" this 
operator 



2. asV ^ 1/s, Ord 1 

3. significantly reduced fluctuations cf. det Q 
Writing the partition function as, 

Z = JvUOrt^ [U] Zla [U] e"'° (9) 

the algorithm will exactly reproduce the full un- 
quenched partition function, if a Metropolis ac- 
cept/reject step is included to incorporate the op- 
erator, Ord 

3.1. Implementation 

The updates of the gauge and boson flelds are 
identical to those used in the original algorithm, 
with the additional constraint that the forward 
and backward rates for updates are identical. 
This is required to ensure detailed balance for the 
entire Markov step. 

The exact algorithm is 

• calculate Ord on the gauge conflguration 

m 

• update the bosons and gauge flelds accord- 
ing to the local boson method 

• recalculate Ord, and accept the new gauge 
conflguration, {U'}, with probability 

To calculate Ord fully, the fermion matrix must 
be diagonalised. This is a computationally inten- 
sive step. In practise, for a suitable choice of the 
parameter, e the dominant contribution to fluctu- 
ations in Ord come from the lowest lying modes 
of the fermion matrix. As a result, the eigenval- 
ues product of (8) can be cut-off at some lower 
bound. 

If all eigenvalues are extracted or a stochastic 
estimate of the determinant is used, the fermion 
mass appearing in the auxiliary boson coupling 
matrix need not be equal to the true fermion 
mass. The mass used in the local boson partition 
function can be tuned to maximise the acceptance 
probability of the Metropolis test. 



3.2. Acceptance rates 

The acceptance probability of the global 
Metropolis step is high enough to allow the use of 
low n polynomials. Table 1 shows how the accep- 
tance probability alters as the number of sweeps 
over the lattice between Metropolis tests is in- 
creased. 
Table 1 

Acceptance probabilities for global Metropolis 
step 16 X 16 lattice, j3 = 3.0 (n = 16, e = 0.05) 



# Sweeps 


Acceptance rate 


per acc/rej 




1 


0.50 ±0.01 


2 


0.42 ±0.01 


4 


0.30 ±0.01 


8 


0.24±0.01 



4. Discussion 

The addition of the accept/reject step provides 
a possible method of reducing the constraint of 
the linear rise in autocorrelation times with the 
number of bosons. It also guarantees exactness 
for any choice of polynomial. This is of impor- 
tance in Schwinger model simulations due to the 
appearance of the topological AZMs. 

Whether the method improves performance for 
4d theories is less certain, as the cost of calculat- 
ing Ord (either exactly or stochastically) may be 
prohibitively expensive. 

Comparisons of performance against HMC are 
currently under investigation. The algorithm 
seems to have better update rates, particularly 
related to the topological tunneling problem. 
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